{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deep Ensembles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import torch\n",
    "import torch.utils.data\n",
    "from torch import nn\n",
    "from torch.nn import functional as F\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "from swag import data, models, utils, losses\n",
    "from swag.posteriors import SWAG\n",
    "\n",
    "import tqdm\n",
    "\n",
    "import os\n",
    "\n",
    "torch.backends.cudnn.benchmark = True\n",
    "torch.manual_seed(1)\n",
    "torch.cuda.manual_seed(1)\n",
    "np.random.seed(1)\n",
    "\n",
    "import hamiltorch\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def featurize(x):\n",
    "    return torch.cat([x[:, None], x[:, None]**2], dim=1)\n",
    "\n",
    "class RegNet(nn.Sequential):\n",
    "    def __init__(self, dimensions, input_dim=1, output_dim=1, apply_var=True):\n",
    "        super(RegNet, self).__init__()\n",
    "        self.dimensions = [input_dim, *dimensions, output_dim]        \n",
    "        for i in range(len(self.dimensions) - 1):\n",
    "            self.add_module('linear%d' % i, torch.nn.Linear(self.dimensions[i], self.dimensions[i + 1]))\n",
    "            if i < len(self.dimensions) - 2:\n",
    "                self.add_module('relu%d' % i, torch.nn.ReLU())\n",
    "\n",
    "        if output_dim == 2:\n",
    "            self.add_module('var_split', SplitDim(correction=apply_var))\n",
    "\n",
    "#     def forward(self, x):\n",
    "#         return super().forward(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb2bae4e1d0>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1dkH8N+TsCWAgklcAJNQS8UNRQIuVFFAxCq7UDGgRRQV2uL7qkhFxaKIVVwoi34oCAJRRFkEF6gLIr6oEKhFQVFEEkFkVREhIWSe948zkUmYfe6de2fm9/185gOZ5d4nM5lnzpzznHNEVUFERO6V5nQAREQUHBM1EZHLMVETEbkcEzURkcsxURMRuVwtOw6anZ2t+fn5dhyaiCgprV27do+q5vi7zZZEnZ+fj+LiYjsOTUSUlESkJNBt7PogInI5JmoiIpdjoiYicjkmaiIil2OiJiJyOSZqpxQVAfn5QFqa+beoyOmIiMilbCnPoxCKioAhQ4CDB83PJSXmZwAoLHQuLiJyJbaonTBq1NEkXeXgQXM9EVENTNROKC2N7HoiSmlM1E7IzY3seiJKaUzUThg7FsjMrH5dZqa5noioBiZqJxQWAlOnAnl5gIj5d+pUDiQSkV+s+nBKYSETMxGFhS1qIiKXY6ImInI5JmoiIpdjoiYicjkmaiIil2OiJiJyOSZqIiKXY6ImInI5JmoiIpdjoiYicjkmaiIil2OiJiJyOSZqShzcZ5JSFFfPo8TAfSYphbFFTe5Us/U8fDj3maSUxRY1uY+/1nMg3GeSUgBb1OQ+/nZpD4T7TFIKYKIm9wm3lcx9JilFhEzUIlJPRFaLyH9FZIOI/D0egVEKC9RKzsriPpOUksJpUZcD6Kiq5wI4D0BXEbnQ3rAopQXapX3CBGDrVsDjMf8ySVOKCJmo1Tjg/bG296K2RkWpjbu0E1UjqqFzroikA1gL4LcAJqvqPX7uMwTAEADIzc1tUxJspJ6IiKoRkbWqWuDvtrAGE1W1UlXPA9AMQDsROdvPfaaqaoGqFuTk5MQWMRER/Sqiqg9V/RHAcgBd7QmHiIhqCqfqI0dEGnn/nwHgCgBf2B0YEREZ4bSoTwGwXETWA1gD4C1Vfc3esIhchItBkcPCqfpYr6qtVbWVqp6tqmPiEZjt+OajcFRNZy8pAVSPLgbFvxeKo9ScmWjFmy+eiZ4fKs7xN52di0FRnIVVnhepgoICLS4utvy4lsnP97/QT16emUgRSs1FgwAzIcOOWt94nouOlZZmPsxrEjETb4gsEqw8LzUTdaxvvlgTfSTieS46Fp9/ipOY66iTTqC1JMJdiS3QokF2LLkZz3Mlq1i6jgJNZ+diUBRHqZmoY33zxZroIxHPcznFzj74WMcjOJ2d3EBVLb+0adNGXW/OHNW8PFUR8++cOZE9NjNT1bz1zUXE/BvpsaI5V2amtedwkt2/X15e9WNXXfLyrDk+kUUAFGuAnJq6iTpWVYneN0nblUhj+VBxOysSabDnp+Zr4/vBSuQiwRJ1ag4mWomDTbGJdWA3VFUMXx9KEBxMtBMH+2ITax98qDpnDgZSEmCijlUqDPbZKdZEGuqDkoOBlASYqH1FU33AFltsYk2k4XxQFhbauzMMZ46S3QJ1XsdyScjBxFiqD5J5sM/tnK6Kcfr8lDTAqo8gfKs3WMblrGg/8MJ9nB0fqCz/I4swUQfirzUUSxkXW9bRC9Ey3bhRddYs1RdeUH3lFdXly1X37bPo+LG8bt7yv31opN/jRJb/UdSYqAMJ1pL21zIK9obmV+DYBHgt9jQ7V4d23qRpOOL35cnPV+3fX3XxYtXDhyM/vmZlRfe6+XwT8wDaBms0A7/oUxiuR5DGFjVFjIk6kECTIfy9YUMl4mBfgdnSDs3Pa/E8Bmpj7NV0VOif8U/diJb6OU7X/9Zrp0tHvKPjxqn266eanW0ekp2tescdqt99F97xw/6ArqnG38IidFdA9Ux8poDqxWmr9IvHXrXrmaIkxUQdSKi+ad+EGqovMlgiYEs7NJ/n1wPoQxilgOqlWKGf4qygifTwYdUlS0zSrlXLPL333qv6ww/+jx/WJVjXRY1Yz8M6bYFNWoF0nZV1hzauX6YZGapTp6p6PPY8XZR8mKgDiaS7ItRU5ECJID098hZbKvK+FpUQHY6nFFAdmF6kh1ErokS6ebPpCqnq1Zg505ssA73WWVmRvz4+fwsL0UMB1VkY8GtM332n2rmzucu110bYl04pi4k6mHC7JUK1qAMlgmhabCnKM3uODqr/kgKqdzScppWzglTkhPigW7dOtX17c9fOnVW//lr9v9bRjC349E37tqZ9Y6qsVH30UdPCb97cfIAQBcNEbYVw3tD+EgHLt8L2wAPmqXngAZ8ug0gSqe/zn5WllSdk6zO4TRvKfq1f97B++WWAE0c6huCNaQF6Hm1NB4hp1SrTaD/lFNVPPw3/uaDUEyxRu2pRpilTgCNHLA8nYh4P8MMPwN69wL59ZsJcnTpA3dKv0O6/U9Fr73Q0zjvOzD4MNcuNW2mFZfZs4IYbgEGDgOnTzXP+q6Iis3ZHaamZcejveff3PHttRR5OxyYM7lSCKW//zpJ4PbOL0Hrw+ThUkY6NuVeh1iNjAr6eGzYAV1wBlJcDS5cCbdtaEgIlmWCLMrmqRR2qpDnel8aNVU87TfW3v1XNzVU94QRzfe3aqtdcY+p6w+p/ZNVHUCtWmOf08stVy8ujPEiIwcJBmK6Z8ovu3WtNzLNnm0O/+GJ49//6a1NK2KCBqQEnqgmJ0qLet8/yUKIiAjRsCNSqVf16VWDtWmDuXGDePODbb819OnYEBg4E/vhHoHZtZ2JOVN9/D5xzDpCVBXz4IdC4cZQHCrRcqtenOBut8CnGjQNGjozyHF7l5UDLlibW4mJz6nBs3w506QJ8/bX5++nePbY4KLkkTIs6kVRWqn78seo995hWN6DarJnq44+r/vST09ElBo9HtWtX1Xr1zMzDmIRRftep3kpt2jTExJgw/POf5pBLl0b+2D17VNu2NcVAs2bFFgclF3Aw0V6Vlaqvv26+ugOmq6S01Omo3G/iRPN8TZpkwcFCLQeQmalL7lyugJmGHq39+1VzcsxrHW2N9P79R/9WJkyIPhZKLkzUcbR8uepxx5n+yC1bnI7GvTZsMC3pq66ycFJIjaoPzcqqNi5QWan6u9+ZFm2053zwQfOu+fjj2EI9dEi1Z089tsqFUhYTdZytWWMGIps1U920yelo3KeiQrV1azPle8eO+J578mTzV//uu5E/9ocfzIdwr17WxFJRoXrTTSaeoUPNNzNKXcESNTcOsEFBAfDee2bQ6corgZ9+cjoid5k0CfjPf4BnnwVOPjm+5x40CDj1VOCuu8LbktHXlCnA/v3A/fdbE0utWsC0acDdd5tjX389cPiwNcem5MJEbZNWnxZhcXovfLv1CIY1XchdP7x27ABGjwa6dgV6947/+TMygHHjgHXrIntJDh4Enn7axN26tXXxiACPPWYuL70EXHMNcOCAdcenJBGoqR3LJdW7PnwHtsbgPgVU59QZxPppVR0wQLVOHQ08SzAOKitVCwpUmzZV/eWX8B5TNfC5YoV9cT33nKkGadtWLav3psQB9lF7xWviiU+pWAXStT1W6nH4Ub9p2t6e8yWIFSvM0zJqlNORqL7/vollzJjQ9z182Ex4uvhi+wf9Xn3VfJAVFKj++KO95yJ3YaJWjXzxHQt2/ai6bEG+NsRPejE+iLmGN1FVVKiec45JeOG2Yu3Wu7dq/fqq27YFv9/MmealfO21+MS1ZIlZzOmii0wpH6UGJmrVyBZHinW3Fj/negHXKaB6110W/k4JZMYM81S8/LLTkRy1ebNqRoZqhw7mg8Sf8nLV009XbdUqviV0r7xiukE6dHDPBxvZi4laNfR60r5iXfEuQKIf2nmTAqoLFlj5i7lfWZl56goKXFQv7P3G9DxuUED17qs3+L1bVd30kiXhH9OqrrUXXjCHuvZalu6lAiZq1ciSbyRJPRA/b9qyMjNQdNxxqbU+8aRJ5ulbtszpSLxqfJDehikKqM4fvqLaff57ypVaC4f1+syFYS99GvW3sADGj1fX9OuTvZioVSN7I9m4hvTWrWYyzLnnqh44EPPhXO+XX1RPPln10kstbk3H0nqt8fqWoY62xcfaUPbrm2+aDQwqMhpqG6zRHOzU3cgKezMBq/9mPB7Vm282h+LaIMmNibpKuG9um3cUf/NN1bQ01T59kv8r7T/+YZ6+lSstPGisr4+fb0wlOFVzsVUB1Va1N+oAzFJAdR6uDS/pWvEtLIDDh1U7djTVIJY+j+QqMSVqAKcCWA5gI4ANAIaHeoxrE3UkbC7le/JJ8+zff7+lh3WVn34y3x6uusriA8faeg3w+LLcFvrcc0d3E++F+eoJN+navJPPvn1mnZLsbO+2YpR0Yk3UpwA43/v/hgC+BHBmsMckRaK2mcejOniweQXmznU6Gns8/rj5/VavtvjAsbZeQ7TIPbl5ugoX6gHUuE+wpGvztzBVM0mocWPVM85gjXUysrTrA8CrAK4Idh8m6vCUl6tecolZRS7ZvtKWl6s2aWK+slvOitZrsG9M0SbdOEyoWr7c1Fh36RK4pJASk2WJGkA+gFIAx/m5bQiAYgDFubm5cf0FE9nu3aZOt1Gj5Nr8tKpuOprF9UOKQ+vVzdunTZtmfuVhw5yOhKxkSaIG0ADAWgC9Q92XLerIbN1qWp9Nmpj/J7rKStUzzzSVLbbVTbs4kcbDXXeZd+/EiU5HQlYJlqjDWj1PRGoDmA+gSFUXRLP4EwWWlwcsW2ZWaLvySmDnTqcjis3rrwMbNwIjRtTYTdxKhYXA1q1mrdKtW1NuR/dHHzV7Lg4fbv52KLmFTNQiIgCmA/hcVZ+0P6TUdPbZwJIlZsPcjh2BXbucjih6jz1mPnz69XM6kuSVnm6WaT3nHPM8b9zodERkp3Ba1O0BDATQUUQ+8V7+YHNc7lJUBOTnm+2m8/NtW1v69783rdGtWxM3Wa9eDXzwAfC//3vsLu5krQYNgMWLzRrb3boBe/Y4HRHZJWSiVtUPVFVUtZWqnue9vBGP4FyhqAgYMgQoKTHDViUl5mebkvVll5lkvWUL0KmTWWg/kUyeDDRsaHZSIfvl5gKvvgps3242YuAOMcmJO7yEMmqU6Tz2dfCgud4mVcn6m2+A9u2Br76y7VSW2rPH7FJyww0mWVN8XHABMGMGsHIlcNttpj1ByYWJOpTS0siut8jllwPLl5ttmS6+GFizxtbTWWLGDLNP5O23Ox1J6unfH3jgAfMajB/vdDRkNSbqUHJzI7veQm3bAv/3f6Z1evnlwKJFtp8yah4P8MwzQIcOwFlnOR2NTeI0VhGtBx80A4v33GP6ril5MFGHMnYskJlZ/brMTHN9HLRoAaxaZZJfr15mY9hId8+Oh2XLTFfN0KFOR2KTOI9VREPEtKjbtDE7mq9f73REZJlABdaxXJJuwosLJlccOqR6001mkkO3bmaRHje5+mqznGl5udOR2MTmRZestH272bg3N1f1+++djobChVgnvKQ8qyZXxPDVuV49YNo0YNIk4M03gVatgHffjS4Mq33zDfDGG8AttwB16jgdjU0cGquIRpMmputj926gZ0+grMzpiChWTNTxYsFXZxFg2DDTFZKZacr37rzT+TfiI48AtWsDt97qbBy2cnCsIhrnnw/Mng189JH5AGUlSGJjoo4XC8v82rYF1q0zpVhPPmlmNb7hUGX75s2mX/S224CmTZ2JIS4cHquIRp8+wEMPAXPmmCnnlLiYqOPF4q/O9eubKou33jIzAK++2nzNjXfN9Zgxprvjb3+L73njrrAQmDrVzI0XMf9Oner6NUZGjTKle/feCyxc6HQ0FC0m6nix+quzt7+7c5c0rD/UAo/+8T946y3gjDOAP/3JtHTt9vnnJoxhw4CTT7b/fI5LwIWgRIDp082kmAEDgE8+cToiigYTdbxY+dW5Rn93ndLNuGfJ7/H1Y/MxfLiZHdiyJdC3L/D22/aV8z34oPkVRoyw5/hkjYwM05o+4QSz4t733zsdEUUsUDlILJekK8+zSqAyv0jL/0KUiu3YYdYrPuEEc/VvfmP2Zly9OvbNdD0es8HBmDHm2KNGxXY8ip9168z+ChdeaMo9yV3AXchdLJrdSsLcM/DQIdWiItXLLjO7ngOm1vn661WfeEL1vfdUd+4MnLw9HrMDzYcfmgXqb7xRtXnzo6e79FL31XPHxAX18nabP9+8doWFNm7qQFEJlqhFbajbKSgo0OLiYsuPm5Ty800XRk15eaYf1KLH7NkDLF1q1rxetQrYtu3obbVqASedBDRqdLSM6+BB4Lvvqq/GduKJwIUXmoHLa64x9bpJo6o7ybcyJzMzIQYMI/Xww8D99wPjxgEjRzodDVURkbWqWuD3NiZqh6Wl+S9yFQncuWxBUtm505T4bd5sllLdsQPYv9+cVsRMsGnSxJTc5eWZaclNm9q4Y4vTovnATFCqZor5Sy+ZvusePZyOiAAmaneLNkEUFZnaq9JSUzkydmzStfziKpoPzAR26JBZQGvjRrPw17nnOh0RBUvUrPpwWrTVIAlYKuZqCTbzMFYZGWY1xkaNTCVIIu4mlEqYqJ3mO5ECMJvhVc1YdNHKbEkvAWcexqpJE7M7zO7dZneY8nKnI6JAmKjdoLDwaKKorDTXuXAZzaSWoDMPY9WmDTBzpun+uP32KNcEcfk63cmAfdRukUKDWeQ+o0eb5QDGjzcLfYUthapl7MbBxESQYoNZ5C4ej9kdZuFCU8L5hz+E+UA2MCzDwcREkGKDWeQuaWnA88+bdc779zfruIQlgdbpTmRM1G4RzmAW+wLJRvXrm8HFjAygWzdg794wHsQGRlwwUbtFqMGsBNizjxJfbq7p/vj2W9MVUlER4gEpWC3jiEBzy2O5cK0PGyTQnn2U+J5/3vx5DR2qoddASYE1UuIBXOsjCXCwkeJsxAjg8ceBZ2r/FbdVTDx6A6s6bMHBxGTAvkCKs3HjgKsz3sFfKp7Ae+hw9IZAW8hxDMU2TNSJgn2BFGfp6cALh3qjBb5CH8zHFjQ/emPNqg6OodiKiTpRpOjMOXLWcXmNsQTdoBB0x2LsR0NzQ81vchZu3kzHYqJOJFyIieJt7Ficlvk9XsG1+AItMQBzUJnR4NhvcqynthUTNREF5v0m1zFvCybgDixBd9zX+aNjGwkcQ7EVEzURBef9JjfUMwlDhgCPLjnr2K5njqHYiomaiMIiAkycCFx6KXDzzcCaNT43cgzFVqyjJqKI7N4NtGtn9tNcsybJ9s50EOuoicgyOTlmTZCffgJ69QLKyvzciTXVlmKiJqKItWoFzJ4NrF5tyqWrfTFnTbXlmKiJKCq9epnNBmbPBp580ucG1lRbLmSiFpHnRGSXiHwWj4CIyKV8uzOys4HsbNz3QBr6Zr6GEXd78Oab3vsFqp0uKWF3SJTCaVHPBNDV5jiIyM1qdmfs3Qvs3QuBYsbBP6IVPkX/aw9j0yYErp0Wqd4dMnCguY5JO6SQiVpV3wewLw6xEJFb+evO8KqPg1ik3VGnbD+6dwd+vPexY2uqRY5d/bHqZ/Zhh2RZH7WIDBGRYhEp3r17t1WHTV0cNSc3CTEVPA+lmO/pjW++Afov7IfKZ/9VvaY6VBkw+7CDsixRq+pUVS1Q1YKcnByrDpuaOGpObhPGVPBL8koxaRKwdCnwt0+vr74uTV5e6HNwXZCAWPXhRhw1J7fxN0Xcl3e6+JAhwNChZsOB2bMjeDzAdUGCYKJ2I65ERm5Tc4p4Vpa5+Jku/vTTwGWXAbfc4jPN3PfxgHmcL64LElTIKeQi8iKAywBkA9gJYLSqTg/2GE4hj1F+vunuqCkvz3yNJHK5PXuAtm3NNPPiYuCUU2rcoajIfEMsLTUt6bFjU35dkGBTyLnWhxtV9VH7dn9wnzpKMOvXAxddZGYxLl8O1KvndETuxrU+Eg1XIqMk0KoVMGsW8NFHwO3ZL0PFpgqmFKiQquV0ABRAYSETMyW8PmVFeKBWCcb8ci9aYyX+WjLRfFsErPn7rvnts6pCyqrjuwS7PojIPvn58JSUojcW4DVcg2W4Ep3wrnXjLUk0nsOuDyJyRmkp0qCYjYFoiS/QD/PMbuZWVTClSIUUEzUR2cdbG90QB/Aqevy6m/nPzc6w9PhhX5+gmKiJyD4+E11OwxbMQz98jjNw40lvwuOx9vi/SsKabCZqIrJPjQqmznmbMb7wEywszsXDD9e4bzTVGylSIcXBRCKKK1XgT38ypXsLFpgNCDh3gBNeiMhlysqADh2AjRuBDz8Ezr4mP2mqN6LFqg8icpV69UxrukEDoEcPYF/Jz/7vmGTVG9FioiYiRzRtapL1tm3AdfUW4gjSj71TklVvRIuJmogcc9FFwJQpwFtll2JkrSeq35iE1RvRYqImIkcNHgwMGwY8cWQ4irL+mtTVG9HiWh9E5LinngI++wy4+eMJaLlmAtq0cToid2GLmogcV7s28PLLwIknAj17Ajt3Oh2RuzBRE9GxHFg6NCcHWLQI2LsXuPZas+kAGUzURFSdg5srt24NTJ8OfPAB8D//E8OBkmyNaiZqIqrO4c2V+/cH7r7bVINMmxbFARz8oLELEzURVeeCpUPHjQO6dDE7mq9aFeGDA33QDBiQsK1rJmoiqs4FS4emzy3CixvPRW7FZvS5dBe+m7Qg/AcH+0BJ0NY1EzURVef00qHerosTtq3HIvTEz5WZ6DO8Kcpnvhje40N9oMSxG8cqTNREVF28lw6tOfA3fPivXRdnYwNm4QZ85LkAw/4iCGsNOX8fNDUl2BoiXD2PiJzjb3lTP+7HGDyM+zF5sum3Duu4o0b5X5EPcOWqfFw9j4jcyd/Anx9/x2hck/E2hg8HVq4M47iFhSYRz5mTFDvAMFETkXPC7IJIy8zAnAn78JvfmMkw27aFefwk2QGGiZqInBNo4C8r65jkevwt/bBoEXDoENC7t9l8ICxVrWuPx/ybYEkaYKImIicFqjCZMMFvcj3jDNObsWYNcHv2y1BJjpmHoTBRE5Fzouia6P5zEUbXGouZv/TFZAxN2NroSLDqg4gSS34+PCWl6IWFeB1X4x10Qge878pKjkiw6oOIkkdpKdKgmI2B+C02oy9eRilOTbja6EgwURNRYvEOQB6Hn7EIPVGGeuiNBTjUrIXDgdmHiZqIEovPAGRLbMIcDMBaFOC2jOeheflJs7SpLyZqIkosNQYgu+etx4MFr2HWlxdiUmm3pFna1BcHE4ko4XnymqNX6dPVBxeBhBpg5GAiESW1tG9LMAs3/Dq4+C2amRuSZICRiZqIEl9uLo7H/l8HF/tgPspQN65raNuJiZqIEp93gLElNmE2BmIN2mFo+lTow4m1+FIgYSVqEekqIptEZLOIjLQ7KCKiiPgMMPaQJbj/+AmYUXkDnv058db18CfkYKKIpAP4EsAVALYBWAOgv6puDPQYDiYSkZM8HqBbN+Df/wbeew9o397piEKLdTCxHYDNqrpFVQ8DmAugh5UBEhFZKS3NVOY1b26WRf3uO6cjik04ibopgG99ft7mva4aERkiIsUiUrx7926r4iMiikqjRsDChcDPP5tkffiw0xFFz7LBRFWdqqoFqlqQk5Nj1WGJiKJ21lnAjBnAhx+arRgTVTiJejuAU31+bua9jojI9fr2Be65B3j2WWD6dKejiU44iXoNgBYi0lxE6gC4DsBie8MiIrLO2LHAFVeYjXFXr3Y6msiFTNSqegTAnwEsA/A5gHmqusHuwIiIrJKeDrz4ItCkCdCnD7Brl9MRRSasPmpVfUNVf6eqp6lqclSQE1FKycoCFiwA9uwB+vUDKiqcjih8nJlIRO5WVGSWLbVg+dLWrYF//QtYsQIYMcKyCG1Xy+kAiIgCKioyy5UePGh+rlq+FIh6N/EBA8zmuE8/DbRtC1x/vUWx2ojLnBKRe+Xnm+RcU4zLl1ZUAJ06AcXFwKpVwHnnRX0oy3CZUyJKTIGWKY1x+dLatYF584DGjYHevYF9+2I6nO2YqInIvQItU2rB8qUnnwzMnw9s3266PyorYz6kbZioici9fPZH/FVmprneAhdeCEycCCxbBowebckhbcFETUTuVWN/ROTlmZ+jHEj055ZbgMGDTe5ftMiyw1qKg4lElPLKyoBLLgE2bTIVIaefHv8YOJhIRBREvXqmv7puXaBXL7PinpswURMRwYxPvvSSaVXfdBNgQ2dD1JioiSgxWDhDMZCOHYF//AN45RVg/HjLDx81zkwkIvezYYZiIHfeaVbYGzkSOP98MzHGaRxMJCL3s2mGYiAHDgAXXGBW2Vu71pKy7ZA4mEhEic2mGYqBNGhgVtorLzfLopaV2XKasDFRE5H72ThDMZDTTwdmzTLrgfzlL7adJixM1ETkfjbPUAykZ0/g3nuBadPMxSlM1ETkfnGYoRjImDFAly7AsGFmMowTOJhIRBTC3r1AmzaAxwOsWwdkZ1t/Dg4mEhHFoGobr127gP7947/SHhM1EVEYzj8feOYZ4O23gfvui++5maiJiMI0aBBw663Ao48CCxfG77xM1EREEZgwAWjXDrjxRrMuSDwwURMRRaBuXbMWSN26ZhuvAwfsPycTNRFRhE49FZg7F/jiC7PpgN0r7TFRExFFoVMn4JFHzCa5EybYey4maiKiKI0YYTYauOsuYOVK+87DRE1EFCURYMYM4LTTgH79gB077DkPEzURUQyOP95Mhtm/H+jbF6iosP4cTNRERDE66yxg+nTgzDPNNHOrcYcXIiILXHedudiBLWoiIpdjoiYicjkmaiIil2OiJiJyOSZqIiKXY6ImInI5JmoiIpdjoiYicjlbNrcVkd0ASiw/sHWyAexxOggHpfLvz989dbn9989T1Rx/N9iSqN1ORIoD7fabClL59+fvnpq/O5DYvz+7PoiIXI6JmojI5VI1UU91OgCHpfLvz989dSXs75+SfdRERIkkVVvUREQJg4maiMjlUipRi0hfEdkgIh4RKahx299EZLOIbJctJ+0AAAKFSURBVBKRK52KMR5E5EER2S4in3gvf3A6JruJSFfva7tZREY6HU+8ichWEfnU+3oXOx2P3UTkORHZJSKf+Vx3goi8JSJfef9t7GSMkUipRA3gMwC9Abzve6WInAngOgBnAegKYIqIpMc/vLh6SlXP817ecDoYO3lfy8kArgJwJoD+3tc81Vzufb0TspY4QjNh3su+RgJ4R1VbAHjH+3NCSKlEraqfq+omPzf1ADBXVctV9RsAmwG0i290ZKN2ADar6hZVPQxgLsxrTklKVd8HsK/G1T0APO/9//MAesY1qBikVKIOoimAb31+3ua9Lpn9WUTWe78iJsxXwCil4utbkwL4t4isFZEhTgfjkJNUdYf3/98DOMnJYCKRdJvbisjbAE72c9MoVX013vE4JdjzAOAZAA/BvHkfAvAEgJviFx054Pequl1ETgTwloh84W11piRVVRFJmNrkpEvUqto5iodtB3Cqz8/NvNclrHCfBxH5F4DXbA7HaUn3+kZKVbd7/90lIgthuoNSLVHvFJFTVHWHiJwCYJfTAYWLXR/GYgDXiUhdEWkOoAWA1Q7HZBvvH2mVXjCDrMlsDYAWItJcROrADBwvdjimuBGR+iLSsOr/ALog+V9zfxYDuNH7/xsBJMw37KRrUQcjIr0ATASQA+B1EflEVa9U1Q0iMg/ARgBHAAxT1UonY7XZYyJyHkzXx1YAtzobjr1U9YiI/BnAMgDpAJ5T1Q0OhxVPJwFYKCKAec+/oKpLnQ3JXiLyIoDLAGSLyDYAowE8CmCeiAyGWYa5n3MRRoZTyImIXI5dH0RELsdETUTkckzUREQux0RNRORyTNRERC7HRE1E5HJM1ERELvf/ZSNx7OSJVeYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "arr = np.load(\"data.npz\") \n",
    "x = torch.from_numpy(arr['x'])\n",
    "f = featurize(x)\n",
    "y = torch.from_numpy(arr['y']) * 10\n",
    "\n",
    "x_ = torch.from_numpy(arr['x_'])\n",
    "f_ = featurize(x_)\n",
    "y_ = torch.from_numpy(arr['y_']) * 10\n",
    "\n",
    "plt.plot(x.data.numpy(), y.data.numpy(), \"ro\")\n",
    "plt.plot(x_.data.numpy(), y_.data.numpy(), \"-b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prior(net, prior_var):\n",
    "    loss = 0.\n",
    "    for p in net.parameters():\n",
    "        loss += torch.norm(p)**2 / (2 * prior_var)\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 50/50 [15:29<00:00, 18.58s/it]\n"
     ]
    }
   ],
   "source": [
    "n_ens = 50\n",
    "n_steps = 20000\n",
    "n_data = len(y)\n",
    "\n",
    "de_all_preds = []\n",
    "for i in tqdm.tqdm(range(n_ens)):\n",
    "\n",
    "    net = RegNet(dimensions=[10, 10, 10], input_dim=2)\n",
    "    optimizer = torch.optim.SGD(net.parameters(), lr=5e-6)\n",
    "    criterion = torch.nn.functional.mse_loss\n",
    "    noise_var = 0.0005\n",
    "    prior_var = 100\n",
    "\n",
    "    for epoch in range(n_steps):\n",
    "        optimizer.zero_grad()\n",
    "        preds = net(f)\n",
    "        loss = criterion(preds, y) / (2 * noise_var)\n",
    "        loss += prior(net, prior_var) / n_data \n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "\n",
    "    preds = net(f_).data.numpy()[None, :]\n",
    "    de_all_preds.append(preds.copy())\n",
    "\n",
    "de_all_preds = np.vstack(de_all_preds)\n",
    "# np.save(\"deep_ensembles_preds\", de_all_preds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1, 4)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3gU5drG79m+6SFUKYmIBRVR4dgVG4pYUFHhiL1gPbbj0WP57BU7ehQQC0q3ix1QwU4RVARBpdeE9GT77vP9cfMym81uSNmQbHh/17VXkt3Z2cnMzj3PPNUQEWg0Go0mdbG09AZoNBqNpmloIddoNJoURwu5RqPRpDhayDUajSbF0UKu0Wg0KY4Wco1Go0lxkibkhmFYDcNYZBjGR8lap0aj0Wh2TDIt8hsBLEvi+jQajUZTD5Ii5IZhdANwKoDxyVifRqPRaOqPLUnreRbAbQAyEy1gGMZIACMBID09vd8+++yTpI/WaDSaXYOFCxduFZEOsc83WcgNwzgNQKGILDQM49hEy4nIOADjAKB///6yYMGCpn60RqPR7FIYhrEm3vPJcK0cCeAMwzBWA5gK4HjDMCYmYb0ajUajqQdNFnIRuUNEuolIAYDhAL4UkQuavGUajUajqRc6j1yj0WhSnGQFOwEAIvI1gK+TuU6NRqPR1I22yDUajSbF0UKu0Wg0KY4Wco1Go0lxtJBrNBpNiqOFXKPRaFIcLeQajUaT4mgh12g0mhRHC7lGo9GkOFrINRqNJsXRQq7RaDQpjhby1sqkSUBBAWCx8OekSS29RRqNppWS1F4rmiQxaRIwciTg8fDvNWv4NwCMGNFy26XRaFol2iJvjdx1lyniCo+Hz2s0Gk0MWshbI2vXNux5jUazS6OFvDXSo0fDntdoNLs0WshbIw8/DKSl1XwuLY3PazQaTQxayFsjI0YA48YB+fmAYfDnuHE60KnRaOKSUlkrfj8fWVktvSU7gREjtHBrNJp6kVIWeTAIVFcDPl9Lb4lGo9G0HlJKyNPSALsdqKgAwuGW3hqNRqNpHaSUkFssQGYmRbyqqqW3RqPRaFoHKSXkAOB0AunprI/RLhaNRqNJQSEHKOR2O1BZCYRCLb01Go1G07KklJCHQrTCrVbTxVJdDYi09JZpNBpNy5FSQh6JMHMlGAQcDlrmXq92sWg0ml2blBJyh4PWuN/Pv91uulg8Hoq7RqPR7IqklJCHwyx0BGiJ22y0ykMh/h2JtOz2aTQaTUuQUkIeCJhulHCYfzudtMx9Pu1i0Wg0uyYpJeRvvQVcey1dKeEwXSwipovF6zXdLhqNRrOrkDJCXlYG3HgjMHUqcPnlpmh7vRRxl4ui7vPpqs9dGj0iT7MLkjJCnpMDTJ/OtMOPPgIuvJCphx4PBd3loptFuV90SuIuiBqRt2YNvwBqRJ4Wc00bJ2WEHAD69mU315wcYPZs4IILgJIS9l6JRCjmdjuFXfvL2zjxLG89Ik+zi5JSQm6zAfvuC7zyCtChA/Dtt8A//wmsW0dBt9mYogjQMtcpiW2URJb3mjXxl9cj8jRtnJQS8uxsoH17oFcv4LXXgK5dgUWLKObLlwOlpbTKHQ7TxaJTEtsgiSxvqzX+8npEnqaNk1JCXlpKt0luLs/NCROAPfcEVqzgDIYFC9h/xe2mv9znYzBU08ZIZGGHw3pEnmaXpMlCbhiGyzCMeYZh/GIYxu+GYdyfjA2LRyRCMQ8GOSWoQwfg1VeBAw7guT1iBDBzJoOgTiddLbqEvw2SyMJWI/H0iDzNLkYyLHI/gONFpC+AAwEMMgzjsCSstxbp6Qx0FhfzTjozk4L+yivA0UcDW7cCl1wCvP22mckC6BL+Nkddw6lHjABWr+ZVf/VqLeKaXYImC7kQNebBvu3RLMl/djvQrh3QqROt7rIyPme3Ay+8AAwZwuevvhoYO5YCnpbGEn49VagNoYdTazQ1MCQJCdeGYVgBLATQC8D/ROT2OMuMBDASAHr06NFvTaIMgzrYuJHZKZ06UZSLixnUFKEBlpYGPPUU3S0AcMMNwJ130nKvqKDg5+YyY02j0WhSDcMwFopI/9jnkyJpIhIWkQMBdANwiGEY+8dZZpyI9BeR/h06dGjU57hcFPD16+kqycujMDscFPKSElZ/3nsvDbXRo4GrrgK2bKHIBwJcRlvmGo2mLZFU21REygB8BWBQMteryMlhyqHLBRQW0o1iGLTQ27Xj8+XlwODBwEsvMXvlgw+As88G/vjDdLMUF2ufuUajaTskI2ulg2EYOdt+dwMYCOCPpq43Hh6PKdzKXVJZyZ/Z2cBuu/FnKMTCoddf53OLFwOnngp88QWzWUSY/eL16lJ+jUaT+iTDIu8C4CvDMH4FMB/ATBH5KAnrrUU4TIu7uprWudNJcS8qovsEoMh37MgMlw4dgPHjgUMOYUbL+ecDTzzBuhHD4AWgulrP/dQ0Ed2oS9PCJCNr5VcROUhEDhCR/UXkgWRsWDwiEYqwmt2Zl0fL3GplBsuGDfSBA/Sdd+3K9MRHHmFCQygEPPoocPrpXNZq5YXA44nTaKupJ+fOPLl3USEJBnm8W7R6Vzfq0rQCUiZ/4/ffgeuuA/76iwIsQgu9UyeW7aelUYyVxV5Rwfe1b0/r/LLLgMceoy/9xx+BAQOAadOofcEg8863W+dNPTl35sm9CwvJihXAhx8y5tFiLjLdqEvTCkhK+mFD6d+/vyxYsKBB7zn/fGDKFArvGWdwwMRee1Gk09Mp3Fu28Byy27lcKARkZNAFU1pKkd+6ldks333H9Z54IjBqFLD77nyPYQD2/feGc+0KGLEbkZ/PIpMdUVAQv4FTfd/fEHbmZ7USAgHgySdZ/+PxAMcey4tyx44tsDEWS/yriGHoRj+apJMo/TAlhHzFCqB375rnhWEAgwZR0A86iP7wUIi55lVVtNp9Pp70GRn0qXs8fK2iAvj8c2a2VFXxQnDHHTRk09OBcHoWDETggg82hGt+aH1Ozp15crcyIYlEzLz+6J9qE2M31Yi6WobDNTfZYuFF2WrlIxIB3n0X+L//43dCLROJAH360DovKGjERqsWuGvXsvxfVYjWh13wQqppOZo1j7y5eeSR2pokAnz6Kf3dF10EfPwxl+nWjf5xm81MSSwrY7qi283X2rcHhg4FJk8GjjySLpW77waOOorrdHTvBAsi8MINL1xmmWp9u+glWq45uvDtzM+KIRKhW8rnMy+SatiHz0d3VShkHjvDoPBGP/x+ukbWr6eOrllD/VuzBli5Eli2DPjhvs/wn5wxyLdvwLBhpoirbQCA334DDj+cPxtEU11TdbUL0Gh2FiKy0x/9+vWT+lJVJdKpk7Lp6n4cdZTIBx+IeDwiGzaIrFrFn7/+KvL11yI//STy998i69aJLFkiMn++yIIFIo8+KtKli7meE/bbKD/YjxIf7FKBDKlEugTcWSITJ9ZvoydOFElLq7lxhsGf+fn1X09jPystLbmfEUMoJFJZKVJRwUdlpUh1tYjPJxII8PVwOPH7w2GR8nKR1atFli8X+fNPHqstW0RKSkQKC0U++0zkpptE+nQvFgPheh1/QCQjg8e13uTnx19Rfn791zFxIpc3jOQfX40mCgALJI6mtnohF6FQjBol0qFD/U7m/v1FpkwRWblSZP16keJiCsZXX1HQf/+dYrF6tcjixRTzb78VuewyEZeL67AYETnb/bH8ht5S3m1fqXh5qng8IpFIPTdandzRIt5cQrsThSQc5vGoqqJo1yXYsYRCIqWlPC7Ll/NnSQnXs3mzyLPPigwcWPu6ZEFI0AAxt9tFZs+O+uC69k/ssYm+8Go0rYxEQp4SPnJFYSHw8stskLV5846X790buPRS4Kyz6GYpKWHaoc9Hv3nnzvTBVlXx7A2FgFWrgIkTgfffp8/WZgPOOw/417+Anj3pqnG7+b560YZ8qJGImaCRllb/njXhsOl28flq9rx5/33gjTeAWbNqtk7o3JnB7KVzN2MrOjd4Wy0WHsdhoUmwXD2yZmZJWprZZKsNHR9N2yelg50Kr5fBS68XePNNZp+sX7/j93XvDlx+OV2fIrwIlJdTmHJyKCyBAH24bjfP+Q0beK5/+infY7Ox1H/kSE4oysxkFWmioTTbaWXByMYiwv0iUn8RF+GxqqykL90weAFdvpyNzaZO5XEAuB/79GGcYsAAXrTvvNN8PR6qx06igi7DAN5sdwPOKh4PF3ywRDflVEKtfOSJhF6jaUW0CSFXAbRIhNab3w9Mnw4891zicY3R5OYCF1/M4Gh6OkXG76d4A2arW6fTTEXctInZLZ99Zgr64MG8MPTuzQtBdjYvBkatfEW0GYsvEOC+Skurx8ULZutgn4/L22wU5/vvZ6qgYp99KNz/+Afjs+np7I/z9NPJqbi1IIyPMBhH4zukwWOKefSFtClZK/Whudev2WVoE0IeCvFhGDXT2YJBCvoTT9TMaEiE1cr88XPOAfbbj8/l5tZsxmWzUdABVof+8QcHWHz6SQThiAUGIjjW+QNGXGDgsJuPQG4urXSHw8xjB9BmLD6V0qkueokQoXhXVvJ3l4vHZ9Qo3kH5/dyvQ4dSwHfbbVvu/rYL4ccfc7n4X0sBYCAvw4fThrpw1FEU/mefBebNS7xNNgQxC8ehL36DG14YAKRHASJ/rIDdzmPdbLSR469pHbQJIa+spBVttZrpa6pvihL3Dz+koC9aVL917r4789EHDGAuevv2/IzSUlrnLheQ/e1HyB7zGCKbN2EF9sYkXID3MBQBUOn371aK80bm4vTTaaGnpfF9TieF3Zic2haZSjHckTUeifAiqC6EbjdjDhdfDPz8M5c54wz+nZVlGsQ2Gy/Qs2ax0Keur+SJ+23CuNJz0X7jr0D37jDuvQdy3jB88w1rChLdmVkQxv34P1yGV5HmMmB74Vlg2DCI8Bipi3bSaSN3ZJrWQZsQcp+Pj9giE+UGUT8NA/jySw6Z+Oab+q3b5QKOOAIYOBA49FCe2Js2AYHZXyHy0svIxlZ0wGa0Rxm8cKMKWXgLw/AqLkMZ2gFg75ezz2YVqmp94nRS0NLSmtnya0aqq/kzPT3xMiqg6fXy/3W5WLxzzTX0c3frRqu8oMCcqRqJ0GVjsfA4Pfhg4rCBxQLccdqvuP3zE+D0l8OKMCyIwIiybn0+ToeaMCHxdvawrMfjV63ESQ8dA4fD/A7ZbNzmpF9020iMRNM6aBNCrpBtFYPqoSoCo/8VZbEvXMhb7w8+qL/PtUsX4Pjjaaln3jACFaU+VCMNDgSRhWp0xBYYMOCGHxmoxBc4BWP3fQFLl/L9BiI4Ct9gRO6nOOHuQ+A892xYrfSlp703CcbdqWOdh8P0CrhciTN1gkGKs/Kh2+3Uwqef5usnnUSRVn52t9t0wTgcvHu65Zb4Az+UtX7zzcBDb++DtHXLay+Ul8co6rZ9Ova4qbhmwmF1WvZuN9CvH+8OTj6Z25b9xXS0v/XS2m6Qiy8GPvlEV35qWpw2IeR+f90DIWKFPRw2LfeVK2mpvfMO3Sb1pQB/41h8if3wC9JRAT8csCOCbJTAAgtsCKNdngN5M9/C8le+wZQXy/CpnIwgHACAPBRjyDElGHr7nuj79zvIvO1auH0lsGPbVaWV+0u9Xu7HjIz4rwcCFPFAgBa7z8c7kpkzeSG96y5g+HBWb4bDXMZmo5VusdD1cvXVfF8sWVkMmPbrB7z1FpDXMwtWhCG1u+Bsx4DAcLvx5j8/wtWv1m8GuMPBz7hh6UgMKX8dbsR8yaKDMkD9jpkKcK5Z07j3azRxaBNCHgw2fExbtLgHAmyaNW0ahUFZ0PWlF5bjSMzBgZiPPJTCgjAiSEfo8pFIG3gM8q4dig4lS+FBJt7HmXgb52ElegHguXyO80Pc5/sPumEjbAhvD7zFWpStxUpXPm/l649FiXgwaG7+6aezQ2VODrVq333ZLx7gcw4HUztDIbqubrqJsY9YuncH1q1jAHnaNKBvXyDr0N6wrl8No47Z3kIph3TPx8u3LMPNN9f//7UgjGdwHUbiNbgQqHvhuizqeAFOJeb5+a3m+GpSjzYh5Mp9Up/0t7rw+ZhLPmcOLfSvvqJ/t74YiOAALMZJ6d/j2Kv2hOOUk1FSAviGDUMWKtAJW9AFm2CBgR9wBKbgfHxkPRvhMJCFMvwXj+IKjIcNYaSjGg7E+HxaicXm95tNx2JTKwMB6pSyxJcuBU47jVk/vXtzOpPqc6NiBQDFW8U6br3VFPloevfmHZTfz8D10KEsEHK/G0cgE2EYiIQimD6dLYy93vr9zwYimI6zcBo+hwv+Otef0Met3SmaZqJNCLnPR+vPbqcwxM3bricivL1fv57tb3/4gW7QH39sWG9rux047jjgzDOB/e85A8Vbg/AiHXYE0Bmb0BWb4OzSGQvHLsC9Z/+KRaEDAADHYTZewHXIQTlyUIY0xPgWWvikFzGzT1yumq8Fg6aIu93cd+ecQzfIUUcxTbO0lPtVBRStVrPzpNsN3HNP/FTRPfdkSuKcOcAJJ7CSV3WvtNtROye7qop+m1i27T9VAPbAAyxCCuzA0AYAJ7z4FoehH36FEesWiVl/XHSAU9NMpHT3Q4W6xVfBtaYUjBgGxaFnT2CPPegSeP55YPZstknt169+6wkGOQv02muB0z1v433L+RAI2qEIheiE3yz9sObKB1BQALzw37W4FaOQiXJ8hRNwMmbiVxyAMuTAhxjfxdq1jf/nkkAoRC2KDXCGQmaFrcvFvO/TT6dADx7MHPA1a7j5Kq1PleOHQhwE8uKL8UW8a1f6y+fO5QXkrrt4AXA6ozJ+RoyggEYi/Pncc3V2HzQMBq+feoqdEe884xccYl+EdCS+BfPDjcH4HKVd+3CDGtrdsAU7Ump2UeI1YGnuR0ObZsUSCrHbXkUFf4ZCTVqdBAIiRUUif/0lsnSpyMKF7KA3Y4bIf/4j0rdv/Zo11Wieh1Xyr/TxMvGar2XGDJFPPxWZNUvk/Ru+kEcz7pVe+EMAERsC8gKulk3oKAFYo1aQ37R/qolUV7MxVjShEPd5URG7F44eLWKxcHOHDRP58Ud2Lfz0U/7+88/sMvnddyLvvsumZZdcEn9/5eZymYMO4t9XXMGmZhs2cFvqpJ5Nw6pfmSLFzi5SiiwpRrbMxWFyBOYIEIm7TQMGbGuS1tCmZC3QkVKza4C20DQrFjWiTd26qyq9xrhclKXp8dDY83ppZVZUcP0lJcxN/+ijhns8DjyQ6YzHHENjbtMm+pJnz2bpvwVhjMK/cQ7eRndsgKWFfeTxUg5Vwyy/n/v3scfovwaYK37llWY5vnLFGAaX37CBqZcff8yCn1gcDrpjCguBf/+bd0pz59KST0ujayVesLWhSH4BfGu3wIowHAgiADv8cGIhDsRJmIUgalcFjR/PdgwNZtIkRO68G8G1m2Dv0QWWRx5q8ZiHJvVpEz7y116j79VuN/3kKm/Z5eIJn5XFJJAOHcyTP54NqJ6P/hkOU8BVnrqqaNy6lUKu/PN//QV8/TUfZWX13/6sLHZSPOkkrrOqiqXlb79NMX8Qd+HSrjPR5fFbWvSk9/l4YUtPN5MtPB4zQ+XBB4GxY/naffcB557L15X7w2Yzc8ZLSvj3t9/yfbEYBvD448zlPvVUxiweeoip2wAvAOnp9e+0WCdxfNchWOGFG9/hMAzG55AYb2N2NlMkc3Mb/nFer+n+q6zk9Kq+fRvQOVOjiaFNCHnv3ux5Uh8sFlp27drRR9q1Kx/dujFOlZ9v+m6BmmX+gQDFrLqa4h4M8qQsKqKFHomY/t+//6b1+NNPFK76kJXF4GifPlz3zz9vE3MLpyFddlltt+zOQoQXGHVxjET4PysRHzOGbYStVlrXRx7J9zmdZm64YXAfVlXx7/nzKfjxuOoq/r/TprGAaPfd2XFSFXRlZtZdUdogevcG1tWOPYRhgRduTMy+HreW31Pr9ZEjeeFqCGVldKO/+y4FXOXJZ2fzfx04MEkXJ80uRZsQ8kceYW6xSovz+ymwXi+FprKSj4qKutufKvLy2PO6Vy8GPHv2pJA4HOaYsug2AFYrn1M9X0IhPp+ezu354Qda6YsX1y85oaCAJf2dO7MCddo0Cvi0aXTDtITlpvarahGr3Cl+PxNGHnmEyz3+OLN1VAC6sJDPp6fzParB2fffs+NhvK/ZkCHAbbdxH55yCsVvzBhW1IZCXG9mZhL7oEydClx3HeCtnb4YcaWj6skxGDrlnO2DuRWGweNz0EHxVxsKmS6+SITB7yuvrPtu7aKLgGee4f5SM0mbkoWl2TVoE0JeXk6RUKJqs9X8Gd1rJRRiVtrmzeb8x9Wr+XPVKj6neojE0rEjBb1XLwp9nz7AAQewoRZgWq2lpTxZvV5ug8XC30tKmD732We88NSFxUJBPPdcCvhXX1HgP/mEF5Zma+aUAHWRcrtN95LXywEQN97I//2BB3hHofqxr17NO5cuXcyB114vB1w/9FD8zzn1VPYbb9+elvjYsWxl+957vBOoruYdlRK6pBFdcWm1csO3FenI+SOwZAlw8MG1M6IOPpj9YJzOmtuj0jRVOutVV/HY1Ye8PCbdHH+8eYenhkyrh2q+ptEAbUTIVUGQEuy6lguH+QiFKEahkFm6r4YRbNpEf/fy5XTZLFvGvxPlGhcUsLHWYYfRpdC7NwVdibqy0NVdQSDAdX722Y7z03v2BK6/nimQq1bx1vuNNyh0O6vZViRC95HdzniD38//Y9EiYNgwXkT/8x/gwgsp4g6HWbizxx78f8vKuL/fe88MhsZy8snA7bfT7eXxUMh8PgaSDz3UvIBkZyduDdCc3HNPfH/+//7H/90weExsNnN49JgxvCDV170WzcEHc78ecgi/P3a7aaGL8GLmdmurXdNGhDwRscId/ROoba2rk8Hvp5CoQsGMDN7KK4FfsQL4/XfmH//+e+3qwHbt6AI56ihOcM/O5mdaraaFptq6rljBxl0//ZTY7ZKXxyDfmDHcpv/+l66HzMydI+aqYKd9e+6/sjJmnAwfTsE+80xzaHx6OvDnn9wnvXrRaly/ngL87rvxs1MAFvncdRff3749/8e33mIu+tix5oQm1WQsthhpZ+D18sIaO07Q6eR3oUsXHufKSn4vLrmkfqMHd8R++zFD5sQT+V1wOk2jwu02Lx5K6DW7Hm1KyEVqinW0la0s9ugverR4q/cqVCaKCi6qdcR2UoxEWIY+fz6t67lza7tNDjiAQjVwIHuMuFwUpfJyCmRhIf3nb46pxNrizLj/m9tNUZs+nZ87aRLdEPWdzNNYRFiJabfzwlFaShfRDTewT3jv3nQxZ2QwSLx8OS9Qu+/OC9qaNczE+eor+tHjfa0GDqS/HKBIb9xIfzjAfdqtm+lWSU/no6UyPD74gBeuWA46CJgxg9+nZcvYliBewy/F3nvz7qNzZ+Duu82+7HWxxx70oQ8ezOwrw6CoZ2bWbFPhcKRua2RN42gTQq5GvSkhVsGlaHdLrAAbhmkhq/ep5SIRcwBCvAwCkcS3siLMWPnqK7pOZs+u2QKkZ0+e5Oeea1aS+/0Uy43Db8SM8iPxIYYigtrq7HbTwv/yS2a4fPwxLxI7ms7TFDweWuA5OaZraOxYVmpmZtJH3rkzLcWVK/n/dOtGES8qojvo55/pkojX2Ozkk/maxcLX8/Loc//4Y2atPPaY6R/2+/l/x+vxsjM5/XS6e2I59FCmkf7f/yVu+2IYbM17551mkLy6mgHOZ5+N/x6Lxbz4A7T8r7mGbi0l5h07cn2BAPeVeo+20HcN2oSQV1SYBTtKZKPdFMrytljM4KNKJ1Ri73CYGRXbhwlEi8W0aYzmrVvHFnz33MMzaQf4fLTSP/qIU4pKSszX+vVjL5IBA/j5Vf0PQxFy8T2OwVjcAA9q59dlZjIGt2QJLbT33+fmNEfKmgitabX/ysv5v6jg5ssvA/37U1iLiij4nTqZbWb//psWuuo5HstJJ9GdorKB1HoGD+YxWLKE60pPp9ipIc1JSztsJJs2caZoRUXD3te9OzBxIi/Gdrt5x6iyrd57D7jiivh3LSr4vWwZ71gAxmQeeogxGsPgRTAz0+zoGZ0Oq2nbJBLylCrR93pZGl5WJlJSIlJayt/Ly0UqK0V8PpaRh8NctqKCr1VXiwSDfL6qis95vSy/rvF4c6JE3GkSAcyHO43Px1u2R75EYPBn1DJ+P0vVL7igZqV2VpbIpZeKfJh3kSzCPvIVDpHncal0wOa4JeJ5eSJdu/L3U04R8XjibHMSHh6PyPr1Ips3i6xcKfLFF9xWQOS220RWrRJZt07k999Fvv6aLQyWLxf57TeRDz8UefppEbc7fun9sceKfPutyC+/iPz6q8gPP/AzBg3i69dcw+NRWMifGzaIbN3K49PShEIikyfX3Yoh9nHxxdxXlZXx11dRwddeeCHxOiwWkeuvF7n6avM4OJ1sF/HXXyJ//imyZQu/ZyL83ldU8LsdDu/UXaTZySBBiX7KCXlJCcVb9Vnx+80vb7SAV1Tw9+jXqqr4fDCY4APy8+OfWbF9TxrQS6OsTOSll0T+8Y+ai59q+Vg+xMkyF4fIZJwpXbA+7kfvtpspko880qjdVieRiEhxMcVn5Ur2R9lrL37eoEF8ft06ivDXX4t8/z1FefFikdmz+b9lZsbfbYcdxj4rS5ZwHYsXmxcDJU6rVlG4S0t5bDds4D5LeIx2Ftv6q1QgQ4a6ZtRLxJ94gvuyomJbj5Y4+F6bLBXd95USZMttac/VKeZ33ikyfrzICSeYz/frx/23cqXIpk08ByIR7q/Kytr9cTRtizYh5D4fv6jR4i1St4CL0BKqrOSjzgZbhhH/rDKMmsvVV/C3oazemTNFRowQcTi4uA1BOQ9T5JN2w+XNkXOkXbv4q+3efdvyNp7EycTvp3j+8Qcfp57Kz+rVi1b36tUiP/3EbZ87l4K+bBkF/fXXJeE29+snMm8e17lpEy3JH36gcCthuukmClFhIY9NYSGXrUsIdwpRF+ogrFKIdrI7/q5TxK87cemO7yK2rbcabqlAuqxDJxmCdxOu02YTefBBNm978klzX6eni4waxX2p9h4MFwQAACAASURBVFcoRDGvqGATOE3bJJGQp1SRsNNJv6nDYfq//X6zfFzlP7tcpi9Z9c42jHpkftS3/WiiFrMJnjcMBiqPPZYFIN9+y1YqEYsN0zEcl9qnoPyAYzBqVPy86XXrGPgKhYB//tOsokwGPh994uEw/boff8x9PH489+H69fQRp6cz0JaWxvesXMkGV9GxAMUBBzAfHuA2V1TQ32uxsLfM7Nlcz+23m7n3gJm6uaM6gXoxaZI5AbuggH/Xl7vu2h5xtCGMDijBzzgQB1l/ibv4YfgOj3wzAK53dvAZ29brgg+AARcCuA93oad1ddzFQyHGHf74g9kvU6Yw5766mmmpzzxjps9WV9NXrgKhml2MeOre3I+mtrGNRGidV1aaFng8C87vN10w9bLw6usyaaBFHktlJa3gr78WOfxw8+0XXCDy/PN0OcRbfU4Ofx5xRD1au9aDYJAW94IFIlOmiNjtXP+LL9IF8tNPbD27eDEt6tWr6Xp57TW2nY23jb16sV3vrFl8/7JltOLnz6cFqXz+jz7K41ZcTN94RYXIxo10s/h8TfzHmtpGNs6dWRiGFCFHnnxSxA7/9pe6YY1sQJ5E6nP8o9brh13KkSHr0VkmYNj2fR/vkZ4u8vLLIh9/zP14++3ma+efz9bLmzebsSJtlbdd0FyuFQDdAXwFYCmA3wHcuKP3NEXIAwHzy+rxJA7uKFdLg4Nm9ek9nYR+04EAA1Zr1ojcfTdvowGRgw8Wue8+EatVEp7UgMjgwU0/WUtLKQyzZ4t06cL1XnKJ6QP/4QcK+qpVIn//TcGvy53SowfX9/PPXMe6dexdvnw513HVVVzuwAPN41haav7ctIk/m+wfb+KFNtH7q7rvI2VlIrNxlPTBIumO1TITx0gQ25qyx7rgdrDearilEO1kSccBcuWVst2dEm/TO3Vi4HXWLO7jUaPMZU86iftc9YkvKoofbNWkPs0p5F0AHLzt90wAKwDsW9d7mhLsrM8wCY+HyzXZsquLhg4biEMkwsDehg0in38u0q0bj8juu4v861+JLTR1DRkxovFZCn6/yIoVInPmMCipBHbRIgbSFi2imK9dS4tv3jyRcePMu4J4fvw//6TVroKbW7fSyl6+XOSTT0yhWrCAQl5Rwf+/qorikzT/eH1jHYlIcKEOvTGJvvzd+spC9JE5OFQKEbVDdnShiFlvGIaUuzrKhtFvyezZIgccwJcSBY/32ktk+nSRb77hPhwzxrywH3UU74CKihg03rSpmb//mhah2YS81gqBDwAMrGuZxgq5CugkIhIxJwep1KxUwOfjrfHChebJ3K4dp+4kCoK5XPz96qsbPiEpEuGJ/v335me0a0dRV4HPxYspyosWUSDGjEksMN2702IvKqIF//ffvNsoLeVFYckSulwApjOKmBlEKh1v0ybuA48nCTu0qRa5SMILdTgsUjl+qqy07ynL0VPCMMyra30u5jHr9b8+WUpKuI9eecV0q+Xlxf8XDj9cZOpUM3towgQzRbF/f4p8URH3/+bNOh2xrbFThBxAAYC1ALLivDYSwAIAC3r06JH0fzBaxFPRPxgO05JavtzM6khPFznxxPgndHq66dM+99yGWV/V1fRb33or32+3i7z9Nk/8P//ka3/9RQGfM4c+81gDdbuPuBst+8pKvmfZMoqySiVcsULktNO47D778NZf5VOr41VZSddLUVGShKeZR62FwyJV46dIVfd9mnRHpqiqovDOmSNy4YXc3K5dTWs79jF0qMgbb5ij9CZONIV///0Zeyks5N1QaWnTRyFqWg/NLuQAMgAsBHD2jpZtarAzlkikHjniKYC6GK1fL3L22Tw6brdppcfzm2Zk8Pdjj6WbYkcEArS0R482/fCjRtEN8tdfFIY//mD+98yZTBFUcznjGbhLl1KIi4pofa9cSREpKaE433wzl83Kos89EqGLrLLSFPLCQi6bjADudpLg+qoLVUyVDMJh7oelSzm3VLnYzj47fqzEYqHr7c03ebx++41WeqdOpph/9RXdWyq1U4t52yCRkCcl/dAwDDuAdwBMEpF3k7HO+qKGH4gwpS2VmwipFMlOndjn5Lzz2Ilv+XL2NYllyxamJXbowIEWRxzBjoSJEGGq4aefMvUvHGafk8suY38Vn4/7b/NmpjiOGsW+IPG6NfbsyVTF3Xar2btGDZqorGQ/8meeYQbgSy+xrzvAlFCbzVyvx8P3JLWXzIgRbJQeifBnkkfnRTdiayoWC/ebGlF45ZV8/pNP2LQslkiELY7XrWN7A4uFTc2efZbfnSVL+L7ly80Sfq83fg8cTdugyUJuGIYB4BUAy0Tk6aZvUv1RX1ARs19zW8BmYzOq8eOpP34/e5NkZ9de9s8/mSbdsye7M/bvz74s8aiq4ms33cT9dsop7OGhBk3b7fycX3+lEMycGX89e+7JnOZ27XgM1FzTqiqKdHk5273ecguXv/de5j+73ebkJTXguaqK68jM3LV7bTsczNXPyWGv+6OP5n5duhS44ILay5eXswdOWRk7T7pcpph37Mj9f9NNzEFXA7O93vpNrtKkHsmwyI8EcCGA4w3DWLztMTgJ662TWEu8RUW8KcUnCVDzKsePZxGQargUrzHS/PlsJ3v88RTks87iSVxcbC7j99PKv/Za/n7WWexsaLVSvCMRWuETJlCAE002Ouww9hXr3ZsWZEaG2Q8+PZ3bV1TEC4HHw1F2w4dTSDwecwCH31+ze58eSMwLXbt23IfXXMOfn3/O/uSHHFJ7+VWreExLSnghzchg4dDo0aZlfsstbMAVCpnHQIt5GySev6W5H031kYfDZsl9i0flmzmwJkLf8ZAhZpAzUXbdWWexCZXyaWdmsl/Hs8+KFBSYy40YwYyHtWtFfvxR5Msv2eTrkEPir1c9rr66dmZJIGBmumzezACbCm727k3f96pVPFZlZfTFBwJch2qAVlycmgHq5qC6mtlL330ncvnl3I97783jtNtu8Y/LCScwvfOXXxis/vlnkbfeEmnfnq8feCCzXKqqWtF5o2kUSOAjT6k2toBpiQO0xFt8EnlBAe9tY8nPp282SZSXs5vu55+b7WPjccEFtN7eeoszJqPp0oXtU087jdaf10sLefVq9tYuKoq/zrQ0WnnnnkurLj3d3O+VlZwiZLHQ9TNuHDv/ZmbSb9+lC9+flcUycquVn11VRRkKBvlcZmYrOJathNWrzX161lmMhTz4IN1mZ59de1IVwO/G9dfTordaeVxWreIM0dJSzkN9+WXOoA2HzXjMruzOSkUStbFNqVOn1Yk40OC+K40lMxN48026NtQ4tnhMnMgZm7ffTlE94QS+58kneRE44wwKqMfDW/IpU+huSSTi++7LC4IScdXnHTD93IEAXSPz5rGVO8D5lqpFTXo6BVuEvmDVT149VH8VDenQgfvTbjfjDE8/zcEe//tf/H01bRrw6qumPzwtjTbGSy/xAjt/Pi/Wq1aZgWYVX9KkPil1+qhmQK1GxIH6N9pqIsrinTyZ/vCtW5kxEo/XXuM4soIC4L772MDqjDMopoZBf+m8ebTWJkxI7DO96CJOKdprL74XoBAr1LBl5V+/9FKu+/rrOV0nFKL/22bjslarOS1HoQaBtBmSEC9JS+M4Pb+fwzcOPJBW9VNPcdLSvffGf99rr9FnHgrRv+50Mgj+9NO8KMyYQatciXk4XPeYOk3qkFKnkMvVykQcMKcRR5OWxueTjMPBjITp0+mq2LiRWhGPsWOZoibCoKQKDP/yC8eP3Xhj4puG7GzgxReBxx/n56hhyNEjxUTMbJdAgCJeXMzA3MMPm1OZ3G7+roJtwaAp5Oq2vq1kG2HSJGDkSLraRPhz5MgGi7lh8DjbbNxf997L56ZMYTbKFVdwSHM8xo6lWIfDDH5arbwQ3H47Xx89mumna9Zw/SrwrElx4jnOm/uR7IKgFqeZi0+iUVWRM2aYxSKq/D3eIy9P5JxzOGGod++6g5mqBHzuXBYNVVSYPVBig5E+H8vx5841q0/32otNwDZuZEXnmjUMaG7ebPZRiX6oCU+pXMRVg2S0BthGMMhWB99/z0pbFew+/ngGj1evZnA73scZhshdd7FK+M8/GTz9+WeR4cPNgPnkyWZFbqpWQ++KoC30I2+1JLP4ZAe35lYrb5MHDODAYoCBsX32ib+64mLg7bdphS1blvhjbTb2F3/pJfq0vV4GWAMB+udj0wMrKhiEe+45YNYsBtlee403I16v6Qqy23knkZ3N9WZkcH2qcMtiaWV3WE0hifESm437KT2dx+Cmm3h38+WXDCJbLMCjjzL4GYsIX3v7bVrm6obxuuuYo15dTd/7ypU8hiJ0seiCodSlrZxCbYN63po7nbzVvvpq5mgrF0fv3o372N69OTD6sst4m+1wmG6s7GxTCCIRikpJCasGX30VeOcdbs+ECUDfvvzdamVhiyrysdvNYR/KnaJEo035yJMcL3G7uf+tVv686CI+P2oUxdjp5MV86NDa741E6JL56CPua5eLbrL//IfHe/NmivnmzbzYh0K6YCiVaSunUNsgajLNdjwePh+FYfDEFGEWQ58+9Jc7ncDBB9f/4zp04Mk8Ywaw//4U3I4dmTLYpQutbLebftrqaj68Xp78b7/NDBmLhT8HDOA2eb3mRCTlE7fba6a5qawVoA35x4Gkx0vcbu47VdF77rk8ZsuWAe++awb/n3uOweVYQiHGQ7780kz97NaNYt61K6tG77yTF+aqKl7EdfAzNdFC3ppowK25zWamqL3zDi3gxYspxDfeCBzYY2vCj8nPp4CPH8+884wMrsftptWmhFeNEFMpbU4nT/Rp0xg0AxhcO+UUMxvF66WgO52m0MS6ZUIhCnmbcqsAdKmNG8cdbBj8OW5co11tVivvjmw2Bqzz8hhUBphOumEDj09VFTBmDCt7YwkEgDvuABYupKvM6QR69aKY5+QAc+cyqF1WxuOsxTw1SbmCoDZNA4uLVD64COeAnnoqrd2LjlyOft89i7Xohr+xJ4qRBxf82PP4Agy4al/ssQfTF0tLede/115cnxJYm82ci2qzmbnjFRWcFTl+PJd/8EHgv/81xbuqiuvs2JEXh6oqswAoGq/XTIV0u1O70Vlz4/HQ9ZGezuPz999M7/z5Z/ZjGTeO+zM9nfvxn/9kamksublMQ+zZk8eqpITrePBBCvf11zMdtWNHHhPV/EzTumgTBUFtnkS35oMHxw2AKiEUoWvjkUf4lunf5aMUGeiFZTgTb+FqPI8r8AKO/vkZiFCQc3JomWVkAL//ztv1wkKz4lJd30MhnujLlwPnn08Rt9koCnffTYvPYqFgV1fz5Hc6axYAxRIdVGtTrpVmQOXh+/081t26sbAnK4uFWpMmUaTDYe7radPMLpPRlJbSjbJxI4Xf6aSv/OabeQxeeIF3dkVFZkqiDn6mDlrIWxPxbs0vvpiRxAQBUNUCNRRiZsO55wI+uPA07sKPGAAvXKhANlwIonvZou1Vg506sbDI7eY6fD4KfF4ehcEwgAULWIRy3HHAfvsBn33GAOaECQy0hsN82O18fyDA647dTiFXAh9NtFslma1gdxrN0CCtLtTxVaKak8PGWHfeyb8fe4wX4kiErpOMDLYX7tmz9ro2bWK4pbiY601PBw46iA26AFblzpplVvn6fLryM1XQrpXWTj3dLV6v2R52SOefMNd/KAAgD0U4FR9gd6xEjxwv2k94Bk4nBdnlMgt0li+nIKxezda4S5fWtMjcbt7K33ors1M6dDBT1tLSKA6Vlaw2dTjMrIpYi9znM4uDLJYk9yBvblRWUXRAOi2tSX7w+qC6RqrYg99P18hdd7GnTo8eDD6rfd2uHS/KJ51ECzyWPn2Y+eJycblwmMHTiRP53Isv0t/erh3vBlLqGLVxErlWtJC3diyW+GaRYdTKFVNibn9rMt6/+jM8GLwNS7F/3NWqgKYKSMZitTKTpX9/4JhjOLQiK4vL5+by5FZiDdDaczjYD8Tv5wUiI6O2xV1VxX8pHE5BP+xOapAWSyTCfG/VsAzgsV67lsHq5cspzpMnc99GIry+rFzJZlolJbXXeeihdMWJ8DsTCDBg+skntPpfeYXHXZX6p9RxasNoIU9VGigeHs82S3raNITvfxDvre+Pb9NOxoY9jsHqQFds2UIxVWXy6el0p3TrBhxwAB89e9KyVv3PMzIo3qEQBbpjR7PsPi2NlnhxMRt5ZWZS4G02WnfRqCZbyvXS4n3kG0oDLqrJpqSEdzOqS2QgQFfKmjXML1+7Fjj8cBZlqSwhwwB++okTh6qqaq/zuOPob1cpon4/3Ss//MD008mT6XqxWFLwWLVRtJCnKg28nVdVeirXPBBgpoPfTwFWWmS1crmSEgqxqrq02yne1dUMfobDtMQ7dDDXmZZmtp+NRLhsIEBrHOB64534ajiGw8Gf8Sz2Vk0LWeQAD39ZGe+E1AWyrIyukfXrma1SVAQMGsTxemrsYTBIv/e118a/+xo0iJlIakxfZSVTE3//nb74d94BunfnsdRtb1senbWSqsQGQPPyeDZfeGHcYJvKZFE5334/ravsbJ7UqvJS5YV37MjVOJ208IqL6VetqKAwd+/O9W3dSqFQxUGRiOnrVgJfV5AT4PKqhWp0lWfKsBMbpMXicnF/hcNmd8P0dO7n/Hy6QrKyGJC+/35a4KrZ2SmnsGgoXs7+Z59xPJzKVOnYkf9Ot2502Vx9Nf3zkYhurtWa0UKeCqheLm++SedocXHCEn5VRl9dTVF1OOjz7NSJQhsI8PmKCq4qLY1BrU6dKNxKvLt25W18RgZfz8nhxSAUMvuHqz4q4TDFXmWxxPOnRiJ8qPapKXmbnuSCn4ZgsXC/KpeY6m6Ync3fe/RgCqHLBUydykItr5cWdiBAX/kTT8Rf94cfslNmIMDle/QwO19++y07J6rc/+gWxPVmJ2f67IqklGul0V+ktsLeewPr4lR/du8BLF++XSwBs+1stBWmWs6q6UAqWKlSBiMRLhMMchm1LuUnF6FVpvLDnU4u7/HQP64mtmdk1N7EQMDMhVYFRHpOZ8OoruZdU3o6911aGvd3aSmv7R4PG2rdeivF/e67OQ0qHDYv1s8+aw7/iOXii1nqn53NC64aUBIMcgbrbbeZQ6LrfTc1dSq7dXmjXIPuNPaWGD68qbskJXE6G1/R3CZ85Mqa3GXJyAAQ73gZQFXV9krMRBN31IQlFbwKhWiBqQwS5T9Xt9CqGEUV96iKTxXoVP54u52b5vMlFmjl4rdaU9Q/3goIhejecjq5j9XIPREex7Iyiv3kycBDD3H/PvMM0K8f39ujB63shx+moMfjiis4vEJZ+l9/TescoMvm8sv5+fVOSezdO7HxUVc7zjZMdF//htImhHyXJwnBtmCQgqss6nDYnAEZPbQjGDQFXblK1LAIu53LlZZSoHNzuYyImR4XjQh9tqriU40i0zQcVazjcpnHEDADyXY775Iee4wi7nCwgKtbN/7eqxdF+Oab6X6JxTDYq+fUU3ncwmHg/fe5rNVKT9I555ifv0NaMNOnLaKDnW2B+gTbduCPVH5tdXejshGAqNRFmBaf1WrmhauTNxikMBsG//b5aBEm6pmi3GEqx1n3Vmk8LlfNaUsKh4PPqYlQd9/NeHggYA5grqrimDevlxb5P/9Ze/0iHA04axa/C8EgC4uGD+ff110HfPGFWbOwQ3bSKMRdHS3kqcSOgm0N6GeuBFplkKjUMq+3Zq9wt9vs56JOXmVkBQLmQGaAf8dLcVPvUQaYFvLG43Ry3weDZgYLYPZ9VyKfm0uLfNAgWuhXXcULrs/Hyt0tWyjm8drfhsPspbNgAX8PBLjcKafw/SNG0H3j8Ww7pnUZDy2Y6dMaaS4HiHattCUa4HoRoT9Vibh6Tp2csV0Jo8VDBVWVmDudpqWoslqiA62VlRQZEb4vnvtFUz9EaF2rfe9y0Z+tLq5VVTVdLoWFFOF585iG+uabzEzy+cw01Jtuoi88FpeLetu5M4OsmZks5X/vPb5+++3AHT2nIfOmy2HxVptvjK1zmDSJ/QTWrqUl/vDDOyXTp7WhYlQq9tQYtI98V6CB/kjVPzz6xFeWd31K6NWFIDrwpXzrImYAVWWpRPvmNY1HpQIqYcjIMPvTq4tpdDB582ZmryxcSFF+/31e86uqWB9QXU1RjndKZmSwlN/t5nHt0IFdF//3Px7jIe4v8IL3EuyGzbBEB+K3GQ8qUL59uMguLOrKSGpQ1k8M2ke+K9BAf6Q6+QOB2u4U1To1URGIsq5jo+/Kt67eX17O5dQXV7tVmo6621HZQlarmUEUCFDoVZ96gOI9Ywb75mzeDJx1FnW0Y0dgzz2Zmvjww+ytE0tVFf3tKkhdVMSGWnffze34wHsSjsL3mIujEEGUOq1di3DYLEqrrgb8r0+BXLlj119bRJ1jakxjstFC3pZohD8y2mpWJ74ScyXyKqslGuWbjSfM0b51ZaGrKUMpWQjUylBBZuXeUs200tLMCtCKCoqwimt07gx88AFTETdtortl8WIuW1DAlsZjxlDYY6mooBGtqni3bmUPlkcfBfawrsIaFOAkzMKDuBtBbDvAPXpsP+aqXUDgnofgif0uxRll2NZQbjBlODUHWsjbEo2oPIye/xkr2KptaihU08IDaF3saHCyOolVcZG2xpOHw8GLogp4qrsjp9OctQqY6aXV1Xz+7beBf/yDwc6hQ1n04/HQ/6186LvvXvvzSktphate8hUVbK7234vW4lTMQBAO3IcHcCS+x1LnQQje/8h2C1TFTNzrViACC/ywww8HAtimaolGHLYRVF/35nQpaiFva6hy/kiEP+vhf1QnWjhc25XicJjl99uzFGAGNetCZVBkZNBi1L7x5KIKS7zemndUatanxcJ9ryziUIipiS+/DBx2GN0kw4ax3a0aENKpE/D662zTEEtJidktUTVL6z5sAM6+LBu3p41GR2zGfByCg8Pz8dDK87d3VVTY8rtCAJQiF3444IeT7pg2nIqoYhZNqeasD1rINQDMQc6qajMam41ioLJaVGHQjixsJfbKcteVnMlFZRxZLHSjRLtSVMuF6CyijAxelPPzGaw88ki6Sc46i9d8t5sX3J49KfZdu9b+zNJSzmndupXfA48HyD//GOz+9A2445nOGDwY8IeseOAB4OCDOXpOXWC89z4GcaXDgMCKMAwI/O7cNpuKqOIW6sLanGgh12wn2sKLnddotZol4eXltNzrEnKVpqjdKc2L3c70Q5XSplwp6oIc29LCZqMbZY89gFdfBY46iqJ85pnMagmFWMZ/0EF8PZ6Yl5cD991HH3thIcW9Vy9anQMG0Grfe2/GMocPB449lp0UQ+cMh3vMM8jpkYMIbLB13w2h0S8iPLxtZq2o7K3YvvzNgU4/1NRAWd0iNUv2o19XzRczMszgWizKDROvgZYmuSg3h0oXDQb58HgozLm5td1a1dvSvktLWeH57bdsfDZhArNXVBvkdesoxhs21P5cm43VowcfzPd26gT8+iuDqT17MkPmhRfowklLAx58kDnrhsHPV5Wo0bUMbYV4qb3JQKcfauqFyjgBagc4AbOneWYmv6zV1bWtdzU7VPvEdw6q/7tyianOiNnZPJ6VlWYPeYXqC9+5M2d1Ksv8wgsZAK2qopDn5lLcd9ut9ueGQpxINGECMH8+Y5Z9+tAaX72aF/kJE4AhQ/hd+ve/2ZCruBhwvjMZkd77wZKZhvA++yH0xuSdsq92BsqlogZn7wy0Ra6Jiyo2UZknyupWllRaGrbnCavpMsoPqC4AuoJz56HaBCv3l8Lr5fMqVqFcMMqKV7f9GzYAl1wCzJnD7BZlmavjvXo1cM01dSeY5OcDJ57I0XNlZXTV+HzAPvtwPfffz4tFQfsqTKs8BXv5f4UdvPoYbjfSXn6uTRQH+XzNN8pQW+SaBqEqNpWgqyG9kYiZiaAabqniH1VxmGi4hKb5ULGI2EC1SlOMrhfw+/m7xWLO6+zShdb1iScyO+Wii4A//mDRUF4e88/ff58/E7FmDScVDTm2DO8MeRm9nrkCliULMW8evz9vvAH07Qus3pqB4/2f4gOcgWqkQSAIe/0I3XlP8+2gnYSKS6j9vrNIipAbhvGqYRiFhmEsScb6NK0DJdTR2SqqOZMierRcKGTeUuqhETsX5V6JDW6q3vTKQlSdL71eMziqcpzz8lh2MGgQfedDhwKff26W2Pfqxb+HDKl7W8okB2/gSlxR+QIWfLIZzjFPYMWHi7F4Ma3y0/EBqpGBy/A63scQlKAdQjAQWLs5pacJtYRLRZEsi/x1AIOStC5NK8JqNbsf1mVpq8kxNtvOidJraqNcJrHeUlUFquIWqsGZ6oOirPisLD5efJH55R4PcOml7HSoyuwdDlrWjz8e328eTQAuzMKpGI1/Y/6samyavQgrVgDX5k7CFRiDCKy4Ac/jOxyOrWiPqpyuCF5xTc0S/gsuYCQ1BQR9Z2apxJIUIReRuQBKkrEuTetDWeYq1zwRyh2jy/BbhkTuFZutZv9yu90MaKteKIrsbL72+OPMMAmFGKR85hkKu+pDf8UVwHffAU8+Cey7b93bJbDgJxyJsbP3xRdfACtP/xcuw+s4G9PhhxvXYSxW2A9AlaSh1Oc0Kz4VxcWtvidLS7lUFDvNR24YxkjDMBYYhrGgSI050TSenXwLarEkTjXUtA6U3ztWyFXPeDWeD+Df6qKrxv2p53Ny+Py//01BNwwObn7sMfMikJPDoOiIEcD339Plcpp7FiyISWGKIgAnZs0C7vv0aHx97mhc324KBuArlKIdrsqYhOJyC4KwoRppCCFGDVtxT5bowp8Wy9QSkaQ8ABQAWFKfZfv16yeaJjBxokhamgi/Q3ykpfF5zS6NzydSUSESidR8/b8t1gAAGp5JREFUPhTi84FAzecDAZGNG0VKSkTCYT4XiYiUl4usWiWyebPIq6+KOBz8mp1+usiGDVxXRYWI3y9SWsrnSsZMkz8d+8r1GC0ZqKjx9Yx9WK0iF10kMnOmSJ8+fK6Pfan8hQLZgjypRLpEYt9kGDtlHzYUj4f7Qu2/5gTAAomjqTprJRW56y5zmrGiFVssmp2Hcq/Eq8xVQc9oVCtc1RhNtRx2u1krEAgAgwczGyUjg+1whw5l06yMDHPClMsFWP95HvKevxcP5z6FNeiBR/Bf5KA07naGw/S133YbC4W6dgV+C/bGNcY4BGFDCFb4EWPetsKeLMEg9130IJW6qNd4vEaghTwVSZTM28a7yGl2jMoXjycYdrvZOiEah8Mcnq2ab9ntZs653Q6ccALwzjssIPrxR2DgQJbdqwKxdu24rPe08xD6azVy3ngBd+RPxTLsg0vSpsNqiT9oedEiivnjj/PCMVMG4h77KARhQxB208XSCsfDRSJ0qdS3Pa3P14BZpw0kWemHUwD8AGBvwzDWG4ZxeTLWq0mAHmirqQPlD49FiU3sXFWbjdakym5RYq560ofDzGbp04cj4fbeG1i2jOI+ZQrw229m5WhGBtdTfeYIBP9cjc6yBePKzsNHH1vi9joHgBUr2CL3sce4La8GL8LoIV9BunaHD25Ijx23Y24JvF6zyKouVPquCoY2R/8hXdmZiqghy9Huldg5iZpdFtXnI15loc/H12PHjVVVcVm73cwxd7u5fGkpl1eFRIEAcM45DHIqLBbmmR98MHDggcARR7D4x+Gg0AUCbJd7333siBiPbliLKzKn4b7K/wAARo9mtemORg62BH4//6fY2baxqCZmQNNmdSoSVXYmLdjZkIcOdiaBiRNF8vMZAMrP14FOzXYiEQbffL7ar6mgp99f83mvV6Sykr8HAlymupoBvOJikU2bGBDduJHrLSwUeeIJkeHDRfbdl8HL2Nhkjx4i118v8sMP3CavV2TlSpHHHhOx2+MHQXtihdxiPCWAiM0mMmkSA6+xwdsGkeRzJRhMvH+jUfuxqor7MRjkPg6FGv/ZSBDs1Ba5RtMGUa2I43WfVK9FW+VqcIjLRas8GKxZpVtaaqY2qkpeFehUrRv++os+7wULGBTduNH8zIEDgTvuAA45hF0RP/sMuPVWfkYs++FX/CNtKV73DIfbzaDo4MGN7JCY5LtXEbP/TF29hFS/FVUgFwjwoWotGjtkIpFFroVco2mDKPdKvFt/JdqxLovohmhqOeUWiEToflFtZ9u1M33t8VwfkQiDolOmsIeLaps7YACLiHJz6W+/8YpKVCOz1vYfjHnoPuQQfPAB1//KK8C55zbCxVJQwArRWPLz2QmsgajMnngtngEzxhCJcFuVq0r93dQ8c900S6PZhVDVnLGBTYC+cKu19muxWS2qoldZ7TYbhTwYNIcJu1zxxdVioZ/8+eepo/fcQ/GeM4fPjxlD63x09t1wwVPr/T/jEFRV0Rfv97NVwJQpNUfa1YskZnj5/Wabg3giHj3b1u02+66rv/XMTo1G02AcjvjphgBFRYmyQmW1RD+nhj7YbGYvFzXiTWW17Ii8PDbL+vtv4PLLuf4nn+RUorxLzsTjxh2wo/YVZ/Zsbucll/A9l18OvHzVfHjyeyNiWOtX0ZykDK9QiP93ojYVqvuncp0Eg+YULdWD6M8/WeqxaFGDPrpeaCHXaNoo8YRZEW2VKwtXlfLHLq/cLarfTiBgDp5oCLm5wPjxdKn07s0sluFjj0PF2Rfj7oynYEXtnMlJkzh56IYbeFG68c1/4I51V6MS6QitWb/jHiwPP1zbud7AnPToroaxqYaqr3sgUNOVooqEvF664484AthrL+CRR/h30okXAW3uh85a0Wh2DiobJV7Wh8pgqaoyX1cZGcFg4vWtWiUybx4zWBpLRQVL9FW2yhlniNx4Y/xMFsMQefFFkady7hcbAmwVgA9kAzqLH3ZmotRFE7NWqqu5D2NL8AMBPl9Zyd9Vqf769SKvvSYyeDAzb9T/kZ7O//mbbxr08TWAzlrRaHY9YrNR4r2uCltUNkVs0DMWn8+s6txjD77P4Wh4JkY4DLz8MnDLLdyGI48EuncHpk6tvazLBbzrGwg/3LgUE1CGXOyOlRiHkTgS38MV8TRLQzeVfRIdNFYWeihk5t4vW8bGYZ98wvmnqiDLagWOOYYJMsOGNX2Grc5a0Wh2UTweimYioYvOTlEBUjVEJBGlpQxiqsEUhtG4zpgiwNKlwKhR3M7dd6dPeUmcETXZRhnulP8DYMFLuBGr0RMGwjjNORNnvjgImZnJ7c4Zzy8eidD3XV5On/+SJRw4XVhovs9i4Xi7f/yDj/btzcKsSIQXrK5dG7dNWsg1ml2UcHjH/T1UEBOguMbLdomlosK0WKN97Q1FhIJ4661cZ8+eFMO//qq9bE/8iadxPbJQjtcwEhNxCQQW5OQAZ53FFMXc3MZtRzRqHq3Ko1+5kq0Ifv+dj9ikl8xMivYRR7C6NSuLFrzDwYuL6mNTXc1snU6dGrddWsg1Gk2rZuFCFv4UFgKHH06x3LCh9nInO77EtMAQZOfn4ZsLx+HGj0/angmSng4cfTRw6KEUzPx8oEMHZs7EG/ig+qAUF/OxaRPvNFaupMW9YgUfsQFgt5vrP+444OSTOctUzatVfcmtVvOi6PWaLq7s7MYPn9BCrtFoWj3ffsuZodXVwCmnAN98wwyZWEaOBP73PzO3ffZsumdmzoy/XtXcSnV6DIfpIvH5arf8jcfuuwP9+wMHHUTL+/DDedFQQq1cUU4nLfjofPtg0JyZGtvjpqFoIddoNK0eEZbkX3YZ3T1nnQV8+GF8sR01ilOMooOsa9cCc+cC8+YBv/wCbNkCFBUBJXUMonS5aLG3a0frvVs3WvI9erDT4157mYFfu51irUbnqTmd6nmgtqtJzbpVRUJNQQu5RqNJCXw+4KGHmOpttwNnnMFe6LFYrcDrrzMbJDYjJxCgyKrh4SpIGf28y1XTglb+foultr9/+/AMK+MNfj/XqdajgsRqoHX0nNToAeaqlW1jBV2X6Gs0muTRjDNjnU6mJJ5zDoVv3jxa5rGEw6z6HDeudum+ap+r0i9VOmVODgON7T+fhIz9C2B3WiD5BfC+NnW7myUSNWNOuWSUS8TrNTN83G6+FgyaRUFWq1nJqUTc4eBFQi1TH1dOQ9FCrtFoGobqKLhmDdVuzZqkTrlXovv447xGrFtHITzuuNrLhsPA9ddz2aqq2i0H0tLMgOb2VgVR2y8iqF5bhKrrbgemTwdgCrjDYYqy308hVr1W3G5a5kqcVSm+el5dWFS2SijEdaqRe8lGC7lGo2kYO2FmrMPB/OvRoyl+n33GIGPfvvGXv/9+4OKLmeWiug0CZuMvtYnhsLn9EQClyEEFsmD4qoCrR8LWdz843p2K9PTa1nb0haG62nSTZGSYMzu9XqZQKh95IGBa+CoNNF7vm6aihVyj0TSMnTAzVuVvH300cO21FMLXXqMhfdBB8d/z3nvMJvnww5rWucViukY8HiC4ZiMCsKEQHVGNNNgQhBMBOBCEe91yuK+/HJGJk1FdbU5FSk/nRcHjMUU9Pd0McPr9TJvctAmorOTfalmfz3SphEKNz7evCx3s1Gg0DSPJPb4ToVwUHg9w/PGsAD32WOCafb/E/160YC6OTfjerl2B00+nO2bPPbnJWVlAWRlQuv/R8G7eCguADFTCiSAcCMCOECIwEIIN6NYD9t8Xb2/jqwKhyuWits/n4/ZVVfHv6OZiyo1iGKYVLsLiocaOrksU7GyGMaAajaZN8/DD8afuJHnKvRoKnfHRVLxQ+C5OwRv4+msX+n23GP/Cl3CjEp/j9Ljv3bCBPc/HjKm5iZ06Ae3d76GDMQ+7yTp0x3p0wSa0RyFyUI4slKMdSuBavxkhL1sRhMPbtiOD26RcJ8pNoqYtZWVRwCMR86Eqag2D61D/U7LRQq7RaBqGGpF21110p/ToQRFvhsHf9renwP+vm3CYtww34yk8hrswPngp+uJ7XIrx6I7VmITL4EUdc9e24fEAq1YBq9AewOA6l821lCLvKJb75+RQpLOzKea5ucw5b9eOz+fmAh078u5BtfaNDmxarTXFXbtWNBrNLoXkF6BqbTEcCMAPGwZiNubhMByOb3A7HoAXbqxBd7yOa/AH9m/RbbXZKOzqkZ5O4c/IYDaLcrvk51P4hwyhm6UhaNeKRqNJOYx1a2GDE0HYkYFqvIzLcSR+xA84Gh/gXByN2chDCYZiGpb3c+HPUC/89lvzZIbsiFCIFaR1VZFGs2ZNw4U8EVrINRpN66VHDzjWrIMHaQjDhj2xEk/iZlyN8XgNI5GLYhyeswJ7nD0AZe5eKHCzeGjpUmDjRrab3bqVj3iTklqSrKzkrUsLuUajab08/DCsI0fC4okgCDscCOB89wxsPXUx7n77QIy23YEOtwMHHgj02NYXPCMDOM41B5HFH8CoKoIlMwfWEeehqu/RCAaZYVJezmVLS5kTXllJS7qsjMHMnUFTh0xEo4Vco9G0XrYFUB133gvf2kIEuvcA7rkfN59/IAp3Y8HQAw8ATzzBgQ2//w5UzvkevslTEUQEYeQgs7IKrtfHIveKMIyjj0VBgTlI2unE9uIfNZAaYGC0pISW/JYtFP7KSl4EKiv5KC83hb+8nLnj9SV64lAy0MFOjUaTEqhug6rjoMUCXH458P77FOMpUzjUYfkBQ7G5JAwgss27HkIavHDltYN/2oeorqaIZmfzfZGI2SzLZjOzTWw2+r0tFrPjoco8UU2vrFbzPaqMv6LCFH31e3W1eRFQAyvGjm34PtDBTo1Gk9KoIhvDoCiGw8Arr/D5t95ik6033wSOKvkWbvRAGIIIrChFO1QhA5bi0u2phJGI2VbWZjMn3qvccBGz14qq3oxupqWqNFVfc4vFrADNzDRzxaPtZPV7p0787GSihVyj0aQMSizdbrPv98svM61vwgRg+HDg8vRncU31Y8hEJUKwIgflKEU2KnL3hGfS58h943+wFW9AdfsuwLXXwjJ4MBwOCnAoRCs6EjF7jOfkxBfm6B4qwaAp7Krplir+iR7aLNL46UB1oV0rGo0mZVG9wX0+Ns56/nnztUH4BAPwNfbEcnS0lMN6+skIzPgCrkgFnAjAghDS7YDt8YdQfdLQ7YWqqjFWVRUvFJmZQOfOtNqB2la2+ltVeYZCprArN4wqDLLZOLyisRa5Hiyh0WjaLCIMTi5axAlD06eE4Q/VNH1zjHLsLUuxB/5CbyzFPvgDOSiGtO+M9A+mbx/dpgZDWK1mhovFwm6MmZlmN8XYUntliUc/1DCJ6IZZnTtrIddoNJq4+Hx0i9jtFOBp0zj5Xg1Qjleok4FK7IclOPCqw3H44RzabLGYk37U7+Xl/Fv5zJWgK8GPDpIqEQfMC0P02Dc1mLkxaCHXaDRtmnCYAu7xUFDT0hiwdDopxivzj8Uvm9rjV/TFfPwDv6IvNqNLjXVkZbFj4skns9Oiy2WOiKuqokCrYRNqepAaJReNCn5GZ8AokddCrtFoNHWgOg76fBRMh4MCHw4D3gnT4bnxv3D4ypGGavjhwjrnnlh40XP4AUfgyy+Bv/8212WzsR/6GWcAAwfywlBYyM/IyOC6lTWenW26S6IHSKjgp5JZw+CyrbKNrWEYgwA8B8AKYLyIPJaM9Wo0Gk1DsFgokpEIrehAwHzNeeF5cDvCCNzzELBuNTK65yH/7lvR5awjMNRJ4V+1igMqPvoI+PFH4Kuv+LBYgGOOAc4+m5OK1NxOu9201pWgqy6JyuqO7nwYDie3EEjRZIvcMAwrgBUABgJYD2A+gH+KyNJE79EWuUajaW78frOdbLQrQw2ssFopxF4vXS92O0VWDUveuhWYORN45x3gyy/NXi0uFzBoEHDaafSpp6VR6AMBrkvNHFV+dOVKUa6VptBsrhXDMA4HcJ+InLzt7zsAQEQeTfQeLeQajaYliRVzn4/Wsio4cjhMaz4tjYHSt95i9eg335jr6dYNGDYMOPdcTiWyWMx0SJWL7nBQ/NVFQgVJG0MiIU/GrIquANZF/b1+23OxGzDSMIwFhmEsKCoqSsLHajQaTeOw2ejTVsU8sb1P/H7Tivd4ODzimmuAOXOA5cuBO+8EuncH1q8HnnqKvvTrrgPmz6do5+TQvWIYpqWuxsKl9PBlERknIv1FpH+HDh121sdqNBpNXKLFPBAwS/ZVVaYa0wZQiMNhvrbXXhyI9OefdLucdBJf+/BD4Mwz6XKZMYPWd04OtuenBwK8aLRWId8AoHvU3922PafRaDStGtVnJVrMRSi2KrNEleF7PDXF3elk8HPGDFri119P4f75Z+Cqq4BDDgHGj6f4d+hg9h9vqp88HskQ8vkA9jQMY3fDMBwAhgP4MAnr1Wg0mmbHbqeYh0I1xTwU4vMqSBkMmo21onE4gL59gYceApYsAZ55BujVi26XO+8E9t2X7QMqKznns7Gph3XRZCEXkRCA6wF8DmAZgOki8ntT16vRaDQ7C9UgS/VJUWLu91PMHQ5a735/bcscoD89PZ1W95VXAr/+CkyaBBx+OKtCn3qKgn7LLcDKlcnffl0QpNFoNNsIBCjWKtskupFWOGwGLR0OBjNjs0+CQb6u/O9eL/PQn3sO+OILLjNsGDB1auO2rzmzVjQajablmDQJKCigD6SggH83EoeDj2CQD9Ugy+Oh1Z2Rwef8flrasXNAo900Ph/FfNAgpi7OmfP/7d1PiFVVHAfw7/fNTG/mhULNlEnl4MKNbWYRQtCiINLamEFhzEIomBa5dGG4KBAhgmgRFUwgubDEjSQVmblxWQZSmklDOeZgSTggzrwZR/21OPd077yu+v7cP573vh8Y5r07vplz5zrfOXPuOb/jpinu2NF2825J9chFJFz79wMTE3HXeXraPQf+2yauVX4jiWvX4sU98/Oudz005IZQKhU35j07625wJse9BwbiYRnABfvKlW7x0NhY/MshS+qRi0i4du2KQ9ybn3fHO1CtukBeXHS9a1+jvF5372s1F+C+fG7jfp3Jnr3/mN8fNI8l+gpyEQnX+fOtHW+BX425uOjGx/0N0Ho93rh5eNj1zi9fjuuveP6XQbLmSx67AwEKchEJ2Zo1rR1v0dCQC9+FBRfiyTD3G1CMjLjAnp11xbOSkr8MGme6ZElBLiLh2rPn/4POtZo7nhEf5r7H7cN8fj7efGJ42B2/ciXehCLt9XmFuYJcRMI1Pg5MTgKjo+7O5Oioe97mjc40pAvjSuXWYU7Gy/Hn5lyY37gRfw7/el+cK2uaRy4i0oRkcPubn8ndiPwenXNzboilWnXTFQcG4tfX6/Fq0XZoHrmISAd8zxxwgVypxNu8+TFzvxXcihVuxsrVq64X7l/va5dnTUEuIt0jw8VBaXx4+954X19cQTEZ5rWa6437474nnxcFuYh0B784aHrapaZfHJRxmPvw9j3xxjAH4t6778EvLbkhlzzGxwEFuYh0i5wWB6VJhrevreLL4SbDvFZzC4NIF/xpBbeyoCAXke6Q4+KgND68fV2VZJ2VxjD3W8hVKvksClKQi0h3yHlxUBpfJdEvxU+Ww03e5PTTD83yGStXkItIdyhgcVCa5FL8paX0Oiv+JmkeuwMBCnIR6RYFLA66Fb8Uf2HB9carVRfmyTorlUpcOTFrKmMrIt1jfLyQ4E4zOOjGxut194dAtbq8nG0eW7x56pGLiGSgcSn/zZvLi2Y1bkKRJQW5iEhGGld/+oqJvoKiimaJiAQgbel+sgJiHouCFOQiIhlrXDDke+r9/fnMXFGQi4jkoL8/nlO+uLh8DD3zr5X9pxQREcDNVDGLN3LOa+aKeuQiIjmqVvPf7k1BLiKSs8HBeOaKbnaKiAQoOT6ex81OjZGLiBTAV0LMg3rkIiKBU5CLiAROQS4iEjgFuYhI4BTkIiKBU5CLiAROQS4iEjgFuYhI4DoKcpIvkTxN8ibJx7NqlIiINK/THvkpAC8COJ5BW0REpA0dLdE3szMAwDyKB4iISFMKq7VCcgLARPT0KsmzRX3tNo0A+KfsRpRE5967evn8Qzj30bSDdwxykt8BeCjlQ7vM7Itmv7qZTQKYbPbfl43kCTPryXF/nXtvnjvQ2+cf8rnfMcjN7JkiGiIiIu3R9EMRkcB1Ov1wC8kLAJ4A8BXJI9k0664QzDBQDnTuvauXzz/Yc6eZld0GERHpgIZWREQCpyAXEQmcgjzhdiUHSL5JcorkWZIby2pjUUi+TXKG5Mno7fmy25Q3kpui6ztFcmfZ7SkSyXMkf46u9Ymy25M3kntJXiJ5KnHsfpJHSf4Wvb+vzDa2QkG+XGrJAZLrAWwF8BiATQA+ItlXfPMK976ZjUVvX5fdmDxF1/NDAM8BWA/glei695Kno2sd5FzqFn0K97OctBPAMTNbB+BY9DwICvIEMztjZmkrTjcDOGBmi2b2B4ApABuKbZ3kbAOAKTP73cyuATgAd92lC5nZcQCXGw5vBrAverwPwAuFNqoDCvLmPAzgz8TzC9Gxbred5E/Rn6HB/JnZpl69xp4B+Jbkj1E5jV60yswuRo//ArCqzMa0orBaK3eLrEoOdIPbfS8AfAxgN9wP+G4A7wF4tbjWScGeNLMZkg8COEry16jX2pPMzEgGMze754K8zZIDMwAeTTx/JDoWtGa/FyQ/AfBlzs0pW1de42aZ2Uz0/hLJQ3BDTb0W5H+TXG1mF0muBnCp7AY1S0MrzTkMYCvJKsm1ANYB+L7kNuUq+o/sbYG7EdzNfgCwjuRakvfA3dw+XHKbCkHyXpIr/GMAz6L7r3eawwC2RY+3AQjmL/Se65HfDsktAD4A8ABcyYGTZrbRzE6TPAjgFwDXAbxhZjfKbGsB3iU5Bje0cg7A6+U2J19mdp3kdgBHAPQB2Gtmp0tuVlFWATgU7SvQD+AzM/um3Cbli+TnAJ4CMBKVGXkLwDsADpJ8DcA0gJfLa2FrtERfRCRwGloREQmcglxEJHAKchGRwCnIRUQCpyAXEQmcglxEJHAKchGRwP0LGM8DQstOH7MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# all_preds = np.vstack([np.load(\"preds_de_{}.npy\".format(i)) for i in range(n_ens)])\n",
    "all_preds = np.load(\"deep_ensembles_preds.npy\")[:]\n",
    "pred_mean = all_preds.mean(axis=0)\n",
    "pred_std = all_preds.std(axis=0)\n",
    "pred_upper = pred_mean + pred_std\n",
    "pred_lower = pred_mean - pred_std\n",
    "\n",
    "plt.plot(x.data.numpy(), y.data.numpy(), \"ro\")\n",
    "# plt.plot(x_.data.numpy(), y_.data.numpy(), \"--k\")\n",
    "plt.plot(x_.data.numpy(), all_preds[:, :, 0].T, \"-b\", alpha=0.05);\n",
    "plt.plot(x_.data.numpy(), pred_mean, \"-b\", lw=5)\n",
    "plt.plot(x_.data.numpy(), pred_lower, \"-b\", lw=2)\n",
    "plt.plot(x_.data.numpy(), pred_upper, \"-b\", lw=2)\n",
    "plt.ylim(-1, 4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py37",
   "language": "python",
   "name": "py37"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
